function [ ] = plotTotalCosts( ctrM, ctrC, period, E_c_ref, ...
    filename, titleString, optfig )
% Plot total costs under the emissions market and under the command and
%   control scenario, in one compliance period
%
%   INPUTS:
%     ctrM - Counterfactual with emissions market
%     ctrC - Counterfactual with command-and-control scenario
%     period - Period in which to plot
%     E_c_ref  - Reference emissions level in command-and-control

%% Setup data

% Emissions market
ctrM.period = period;
t        = ctrM.idx_agg_t;
E_market = ctrM.output_tk.Et(t);
Z_market = ctrM.output_tk.Zt(t);

% Command-and-control
ctrC.period = period;
t        = ctrC.idx_agg_t;
E_commnd = ctrC.output_tk.Et(t);
Z_commnd = ctrC.output_tk.Zt(t);

% Reference levels
E_c_ref2 = 170e3;
Z_c_ref = interp1(E_commnd,Z_commnd,E_c_ref); % Cost under command (240)
Z_c_ref2 = interp1(E_commnd,Z_commnd,E_c_ref2); % Cost under command (170)

% Neil Himwich adds unique here, as step function mod doesn't necessarily
% return unique values, but interp1 requires unique values. This won't
% effect normal isoelastic approach, as that will be unique by
% construction.
if length(unique(E_market))>1
    Z_m_ref = interp1(unique(E_market,'stable'),unique(Z_market,'stable'),E_c_ref); % Cost under market (240)
    Z_m_ref2 = interp1(unique(E_market,'stable'),unique(Z_market,'stable'),E_c_ref2); % Cost under market (170)
    E_m_ref = interp1(unique(Z_market,'stable'),unique(E_market,'stable'),Z_c_ref); % Market emissions at same cost
else
    Z_m_ref = Z_market;
    Z_m_ref2 = Z_market;
    E_m_ref = E_market;
end
%% Plot

% Lines
f0 = figure;
l1 = plot( E_market/1e3, Z_market/1e6, 'bo-', 'LineWidth',optfig.axisweight-0.5 ); hold on;
l2 = plot( E_commnd/1e3, Z_commnd/1e6, 'ko:', 'LineWidth',optfig.axisweight-0.5 );  

% Figure formatting
set(gcf, 'Color'       , 'w', ...
    'position',[0 0 500 500]);
set(gca, optfig.labProp{:}, ...
         'Box'         , 'off'             , ...
         'LineWidth'   , optfig.axisweight);

% Axis labels
xlabel('Emissions (metric tons)',optfig.textProp{:});
ylabel('Total variable abatement cost (INR m)',optfig.textProp{:});

% Title
% title(titleString,optfig.textProp{:},'Interpreter','latex');

% Additional lines to indicate comparisons
%   * Same emissions, cost under market
%   * Same cost, emissions under market
plot([ E_m_ref E_c_ref ]/1e3,[ Z_c_ref Z_c_ref ]/1e6,'k--','LineWidth',optfig.axisweight-0.5);
plot([ E_c_ref E_c_ref ]/1e3,[ Z_m_ref Z_c_ref ]/1e6,'k--','LineWidth',optfig.axisweight-0.5);
plot([ E_c_ref E_c_ref ]/1e3,[ Z_m_ref 0 ]/1e6,'k--','LineWidth',optfig.axisweight-0.5);
plot([ E_c_ref2 E_c_ref2 ]/1e3,[ Z_c_ref2 0 ]/1e6,'k--','LineWidth',optfig.axisweight-0.5);

plot([ E_c_ref ]/1e3,[ Z_c_ref ]/1e6,'o',...
    'LineWidth',optfig.axisweight-0.5,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor','k',...
    'MarkerSize',11)
text(1.0*E_c_ref/1e3,1.07*Z_c_ref/1e6,'Control',optfig.legProp{:});

plot([ E_c_ref2 ]/1e3,[ Z_c_ref2 ]/1e6,'o',...
    'LineWidth',optfig.axisweight-0.5,...
    'MarkerEdgeColor','k',...
    'MarkerFaceColor','k',...
    'MarkerSize',11)

plot([ E_c_ref ]/1e3,[ Z_m_ref ]/1e6,'o',...
    'LineWidth',optfig.axisweight-0.5,...
    'MarkerEdgeColor','b',...
    'MarkerFaceColor','b',...
    'MarkerSize',11)

plot([ E_c_ref2 ]/1e3,[ Z_m_ref2 ]/1e6,'o',...
    'LineWidth',optfig.axisweight-0.5,...
    'MarkerEdgeColor','b',...
    'MarkerFaceColor','b',...
    'MarkerSize',11)
text(0.66*E_c_ref2/1e3,0.95*Z_m_ref2/1e6,'Treatment',optfig.legProp{:});


% Axis limits must include zero
% if period == 1
%     ylim([0 25]);
% else
%     ylim([0 14]);
% end
xl = xlim();
xlim([0 xl(2)]);
yl = ylim();
ylim([0 yl(2)]);

% if period == 1
%     text(245000/1e3,5*25/14,'Emissions =',optfig.legProp{:});
%     text(245000/1e3,4.4*25/14,'240 tons',optfig.legProp{:});
%     text(110000/1e3,5*25/14,'Emissions =',optfig.legProp{:});
%     text(110000/1e3,4.4*25/14,'170 tons',optfig.legProp{:});
% else
%     text(245000/1e3,5,'Emissions =',optfig.legProp{:});
%     text(245000/1e3,4.4,'240 tons',optfig.legProp{:});
%     text(110000/1e3,5,'Emissions =',optfig.legProp{:});
%     text(110000/1e3,4.4,'170 tons',optfig.legProp{:});
% end
text(245000/1e3,5*yl(2)/14,'Emissions =',optfig.legProp{:});
text(245000/1e3,4.4*yl(2)/14,'240 tons',optfig.legProp{:});
text(110000/1e3,5*yl(2)/14,'Emissions =',optfig.legProp{:});
text(110000/1e3,4.4*yl(2)/14,'170 tons',optfig.legProp{:});

% Legend
legend([l2, l1] ,{'Command and control','Emissions market'},...
    'Location','southoutside',...
    'Orientation','vertical',...
     optfig.legProp{:});
legend boxoff;


hold off;


% xlim([xmin xmax]);
% ylim([ymin ymax]);
% xstep = (xmax-xmin)/10;
% ystep = (ymax - ymin)/4;


%% Export figure to file

% Crop whitespace 
fig = gcf;
fig.PaperPositionMode = 'auto';
fig.PaperSize = [7.4 7.4];
fig.PaperPosition = [0.25 0.25 6.95 6.95];

ax = gca;
ax.Position = [0.10 0.30 0.90 0.70];

% Save all figures
%[filepath, ~, ~] = fileparts(filename);
%if not(isfolder(filepath))
%    mkdir(filepath);
%end
%fprintf(1,'Writing %s to file ...\n',filename);
%print(f0,'-dpdf','-painters','-noui','-r600', filename);

% Save figure for specific model type and subtype and period
% Determine model type
if contains(filename, "Iso")
    type = "Iso";
    if contains(filename, "FullNoHet")
        subtype = "FullNoHet";
    else 
        subtype = "FullHet";
    end
else
    type = "Step";
    if contains(filename, "NoFreeCyc")
        subtype = "NoFreeCyc";
    else 
        subtype = "FreeCyc";
    end
end

% Get path
[~, old_file, old_ext] = fileparts(filename);
new_filepath = erase(filename, strcat("model/", type, "/", subtype, "/", old_file, old_ext));

if type == "Iso" && subtype == "FullNoHet" && period == 8
    filename = strcat(new_filepath, "Figure_10");
    fprintf(1,'Writing %s to file ...\n',filename);
    print(f0,'-dpdf','-painters','-noui','-r600', filename);
end

end